Optimización Heuristica

Optimización Descenso por Gradiente condición Inicial Aleatoria

Función de Derivadas parciales

partial_dev <- function(x,i,fun,h=0.01){

    e <- x*0 # crea un vector de ceros de la misma longitud de x

    e[i] <- h

    y <- (fun(x+e)-fun(x-e))/(2*h)

  return(y)

}

Para evaluar el gradiente de una función f en x hay que obtener cada una de las derivadas parciales de f en x. Con el código anterior y la función mapply() se puede lograr esto fácilmente, como se muestra a continuación:

num_grad <- function(x,fun,h=0.01){

  # x: punto del espacio donde se debe evaluar el gradiente

  # fun: función para la que se desea calcular el gradiente en x

  # h: es el tamaño de ventana para el cálculo de la derivada numérica

  d <- length(x)

  y <- mapply(FUN=partial_dev,i=1:d,MoreArgs=list(x=x,h=h,fun=fun))

  return(y)

}

Función de Calculo de \(\frac{\partial}{\partial x_{t}} \nabla f^{T}\)

deriv_grad <- function(x,fun,i=1,h=0.01){

  # x: punto en el que se evalúa el gradiente

  # fun: función para la cual se calcula la derivada del gradiente respecto a la íesima componente

  # i: i-ésima componente del vector x con respecto a la que se deriva

    e <- x*0 # crea un vector de ceros de la misma longitud de x

    e[i] <- h

    y <- (num_grad(x+e,fun=fun,h=h)-num_grad(x-e,fun=fun,h=h))/(2*h)

    return(y)

}

Función para Calcular la matriz Hessiana

matriz_hessiana <- function(x,fun,h=0.01){

  # x: punto en el que se evalúa la matriz hessiana

  # fun: función a la que se le calcula la matriz hessiana en x

  # h: es el tamaño de ventana para el cálculo de la derivada numérica

  d <- length(x)

  y <- mapply(FUN=deriv_grad,i=1:d,MoreArgs=list(x=x,h=h,fun=fun),SIMPLIFY = TRUE)

  return(y)

}
optimizador_mult_numdev <- function(x0,fun,max_eval=100,h=0.01,eta=0.01){

  x <- matrix(NA,ncol =length(x0), nrow = max_eval)

  x[1,] <- x0

  for (i in 2:max_eval){

    num_grad_fun <- num_grad(x[i-1,],fun,h)

    H <- matriz_hessiana(x[i-1,],fun,h)

    cambio <- - eta*solve(H)%*%num_grad_fun

    x[i,] <- x[i-1,] + cambio

    cambio_opt <- sqrt(sum((x[i-1,]-x[i,])^2))

    if (cambio_opt<0.00001){

      break

    }

  }

  return(x[1:i,])

}
  • Función de Rosenbrock:\(f(x_{1},x_{2})=100(x_{2}-x^{2}_{1})^{2} +(1-x_{1})^{2},\hspace{1mm}x_{i}\hspace{1mm}\epsilon \hspace{1mm}\left [ -2.048,2.048 \right ]\hspace{1mm},\hspace{1mm}i=1,2\). Alcanza su valor mínimo en \(x_{1}=1\) y \(x_{2}=1\).
f_rosenbrock2d <- function(x1,x2){

  # Esta versión de la función se usa para graficar.

  y <- 100*((x2-x1^2)^2)+(1-x1)^2

  return(y)

}



f_rosenbrock2d_vec <- function(x){

  # Versión vectorizada. Es la que se utiliza para optimizar.

  x1 <- x[1]

  x2 <- x[2]

  y <- 100*((x2-x1^2)^2)+(1-x1)^2

  return(y)

}

Proceso de Optimización con condición Inicial Aleatoria \((x_{1},x_{2})^T\) con \(x_{i} \sim U_{i} (-2.048,2.048),\hspace{1mm}i=1,2\)

Solución

set.seed(123456)

x1_ros_inic <- runif(1, min =-2.048, max = 2.048)

set.seed(789)

x2_ros_inic <- runif(1, min =-2.048, max = 2.048)

x0_rosenbrock <- c(x1_ros_inic,x2_ros_inic)

sol_ros <- optimizador_mult_numdev(f_rosenbrock2d_vec,x0=x0_rosenbrock,eta=1.3)

Graficamente la evolución de la solución:

Función para Generar las imagenes en cada iteración

library(stringr)

library(magick)

plot_generate_function <- function(out_dir,sol_func,name_function,x1,x2,z){

  n_length <- 150

  X <- expand.grid(x1,x2)

  Z <- matrix(z,ncol=n_length,nrow = n_length)

  tamaño_iter <- dim(sol_func)

  n_iter <- tamaño_iter[1]

  for(i in 2:n_iter){

  png (paste0(out_dir,"/",str_pad(i, 3, pad = "0"), ".png"))

  contour(x1, x2, Z, las=1,

        xlab = expression(x[1]), ylab = expression(x[2]),

        main = paste("Función de ",name_function ," iter:",i),

        sub = "Curvas de nivel de la función")

  lines(sol_func[1:i,], type="b",cex=1.5,col="red")

  dev.off()

  #file.show(paste0(out_dir,"/",i,".png"))

  }

  png (paste0(out_dir,"/",str_pad(n_iter+1, 3, pad = "0"), ".png"))

  contour(x1, x2, Z, las=1,

        xlab = expression(x[1]), ylab = expression(x[2]),

        main = paste0("Función de ",name_function),

        sub = "Curvas de nivel de la función")

  lines(sol_func[1:n_iter,], type="b",cex=1.5,col="red")

  points(sol_func[n_iter,1],sol_func[n_iter,2],pch=8,cex=2, type = "p")

  

}

Generar las Imagenes

plot_generate_function("~/Tarea 07/img_ros_grad",sol_func=sol_ros,"Rosenbrock",x1=x1_ros,x2=x2_ros,z=z_ros)

Función para crear el Gif Animado

save_to_gif <- function(dir,out_dir,out_name){

  imgs <- list.files(dir, full.names = TRUE)

  img_list <- lapply(imgs, image_read)

  

  ## join the images together

  img_joined <- image_join(img_list)

  

  ## animate at 2 frames per second

  img_animated <- image_animate(img_joined, fps = 2)

  

  ## view animated image

  img_animated

  

  ## save to disk

  image_write(image = img_animated,

              path = paste0(out_dir,"/",out_name,".gif"))

}

Guardar el Gif Animado

save_to_gif("~/Tarea 07/img_ros_grad","~/Tarea 07","Rosenbrock_Gradiente")

Mostrar Gif


Se observa entonces que la otpitmización por descenso por gradiente con condición inicial aleatoria, alcanza su valor mínimo en las coordenas (1,1), es decir el algoritmo alcanza el mínimo global de la función.


  • Función de Rastrigin\(f(x_{1},x_{2})=20 + \sum_{i=1}^{2}x_{i}-10 \cos(2\pi x_{i}),\hspace{1mm}x_{i}\hspace{1mm}\epsilon \hspace{1mm}\left [ -5.12,5.12 \right ]\hspace{1mm},\hspace{1mm}i=1,2\). Alcanza su valor mínimo en \(x_{1}=0\) y \(x_{2}=0\).


Graficamente la Solución


Observemos el comportamiento de la optimización más de cerca:


Generar las Imagenes

plot_generate_function("~/Tarea 07/img_ras_grad",sol_func=sol_ras,"Rastrigin",x1=x1,x2=x2,z=z)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ras_grad","~/Tarea 07","Ranstrigin_Gradiente")

Mostrar Gif


A pesar de que estuvo muy cerca de encontrar el mínimo global en las coordenadas (0,0), el método quedo atrapado en un mínimo local, encontrando como optimo las coordenadas (0,2).


Optimización Con Algoritmos Evolutivos

  • Función de Rosenbrock
f_rosenbrock2d_vec_inv <- function(x){

  # Versión vectorizada. Es la que se utiliza para optimizar.

  x1 <- x[1]

  x2 <- x[2]

  y <- 100*((x2-x1^2)^2)+(1-x1)^2

  return(-y)

}

Graficamente la Solución


Crear las imagenes para el Gif Animado

plot_generate_function_ga("~/Tarea 07/img_ros_ga",sol_func=opt_ros_ga,"Rosenbrock",x1_ros, x2_ros,Z_ros)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ros_ga","~/Tarea 07","Rosenbrock_ga")

Mostrar Gif


El algoritmo evolutivo encuentra perfectamente la solución en las coordenadas (1,1).


  • Función de Rastrigin
f_rastrigin2d_vec_inv <- function(x){

  # Versión vectorizada. Es la que se utiliza para optimizar.

  x1 <- x[1]

  x2 <- x[2]

  y <- 20+x1^2-10*cos(2*pi*x1)+x2^2-10*cos(2*pi*x2)

  return(-y)

}

Gráficamente la Solución

contour(x1_ras, x2_ras, Z_ras, las=1,

        xlab = expression(x[1]), ylab = expression(x[2]),

        main = "Función de Rastrigin",

        sub = "Curvas de nivel de la función")

lines(opt_ras_ga@population, type="b",cex=1.5,col="red")


Crear las imagenes para el Gif Animado

plot_generate_function_ga("~/Tarea 07/img_ras_ga",sol_func=opt_ras_ga,"Rastrigin",x1_ras, x2_ras,Z_ras)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ras_ga","~/Tarea 07","Ranstrigin_ga")

Mostrar Gif


Se observa que el algoritmo genético es más robusto a mínimo locales y alcanza su mínimo global con satisfacción.


Optimización con PSO


#Función de optimización PSO

opt_pas_pso_ras <- psoptim(

  par = c(NA,NA),

  fn = f_rastrigin2d_vec,

  lower = c(-5.12,-5.12),

  upper=c(5.12,5.12),

  control = list(trace.stats=TRUE, trace=1, maxit = 150)

)
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=14.48, vectorize=FALSE, hybrid=off
## It 10: fitness=2.186
## It 20: fitness=1.351
## It 30: fitness=0.5029
## It 40: fitness=0.2419
## It 50: fitness=0.07325
## It 60: fitness=0.07254
## It 70: fitness=0.008249
## It 80: fitness=0.0001511
## It 90: fitness=1.973e-05
## It 100: fitness=1.973e-05
## It 110: fitness=2.112e-06
## It 120: fitness=4.153e-07
## It 130: fitness=4.153e-07
## It 140: fitness=6.349e-08
## It 150: fitness=1.037e-08
## Maximal number of iterations reached

Valor Optimo:

opt_pas_pso_ras$par
## [1] 1.951811e-06 6.960436e-06

Generar las imagenes

plot_generate_pso("~/Tarea 07/img_ras_pso",sol_func=opt_pas_pso_ras$par,"Rastrigin",x1_ras, x2_ras,Z_ras,opt_pas_pso_ras)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ras_pso","~/Tarea 07","Ranstrigin_pso")

Mostrar Gif


El algoritmo PSO muestra como las posibles soluciones candidatas convergen al mínimo global, siendo este también un buen método de optimización bastante robusto a mínimo locales.


  • Función de Rosenbrock
#Función de optimización PSO

opt_pas_pso_ros <- psoptim(

  par = c(NA,NA),

  fn = f_rosenbrock2d_vec,

  lower = c(-5.12,-5.12),

  upper=c(5.12,5.12),

  control = list(trace.stats=TRUE, trace=1, maxit = 150)

)
## S=12, K=3, p=0.2297, w0=0.7213, w1=0.7213, c.p=1.193, c.g=1.193
## v.max=NA, d=14.48, vectorize=FALSE, hybrid=off
## It 10: fitness=0.3783
## It 20: fitness=0.02858
## It 30: fitness=0.003881
## It 40: fitness=2.138e-05
## It 50: fitness=2.138e-05
## It 60: fitness=2.138e-05
## It 70: fitness=4.776e-06
## It 80: fitness=2.491e-06
## It 90: fitness=2.491e-06
## It 100: fitness=1.403e-06
## It 110: fitness=1.281e-06
## It 120: fitness=1.278e-06
## It 130: fitness=1.278e-06
## It 140: fitness=1.278e-06
## It 150: fitness=1.278e-06
## Maximal number of iterations reached

Valor Optimo:

opt_pas_pso_ros$par
## [1] 1.001130 1.002262

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ros_pso","~/Tarea 07","Rosenbrock_pso")

Mostrar Gif


Igualmente el algoritmo PSO muestra como las posibles soluciones candidatas convergen al mínimo global, siendo este también un buen método de optimización bastante robusto a mínimo locales.


Optimización por Evolución Diferencial


  • Función de Rosenbrock

Gráficamente la solución

contour(x1_ros,x2_ros,Z_ros)

lines(opt_ros01$member$pop,type="p",pch=2,col="red",lwd=3)

Generar las Imagenes para el Gif

plot_generate_function_ED("~/Tarea 07/img_ros_ed",sol_func=opt_ros01 ,"Rosenbrock",x1_ros, x2_ros,Z_ros)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ros_ed","~/Tarea 07","Rosenbrock_ed")

Gif Animado:


El método de evolución diferencial converge igualmente bien, encontrando el mínimo global con satisfacción.


  • Función de Rastrigin

Gráficamente la solución

contour(x1_ras, x2_ras, Z_ras)

lines(opt_ras01$member$pop[], type = "p", pch = 2, col ="red", lwd = 3)

plot_generate_function_ED("~/Tarea 07/img_ras_ed",sol_func=opt_ras01 ,"Rastrigin",x1_ras, x2_ras,Z_ras)

Generar Gif Animado

save_to_gif("~/Tarea 07/img_ras_ed","~/Tarea 07","Rastrigin_ed")

Gif Animado:


El Algoritmo cae en un mínimo local muy cerca del global en las coordenadas (1,0)


Conclusiones